TRANSVERSE INSTABILITY OF PLANE WAVE SOLITON SOLUTIONS 
OF THE NOVIKOV-VESELOV EQUATION 

R. CROKE, J. L. MUELLER AND A. STAHEL 

Abstract. The Novikov-Veselov (NV) equation is a dispersive (2+l)-diniensional non- 
linear evolution equation that generalizes the (l-l-l)-dimensional Korteweg-deVries (KdV) 
equation. This paper considers the stability of plane wave soliton solutions of the NV equa- 
f ^ . tion to transverse perturbations. To investigate the behavior of the perturbations, a hybrid 

^S| ' semi-implicit/spectral numerical scheme was developed, applicable to other nonlinear PDE 

^ . systems. Numerical simulations of the evolution of transversely perturbed plane wave solu- 

Mh' tions are presented. In particular, it is established that plane wave soliton solutions are not 

■^^ . stable for transverse perturbations. 
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1. Introduction 



r^ ', The Novikov-Veselov (NV) equation for u{z,t) = u{x,y,t) was introduced in the periodic 

c^ I setting by Novikov and Veselov [20j in the form 



Ut = —d^u — d u + 3dz{uf) + 3dz{uf), where d^f = d^u, 

with dz = ^iydx + idy) , dz = ^{d^ —idy), where it was derived algebraically from a Lax triple, 
and from this point of view is considered the most general derivation of the KdV equation 

ob'; &j. 

'nI" . While there are no known physical applications of the NV equation, it is related to two 

other (2-|-l)-dimensional integrable systems which have been more widely studied. The 
Davey-Stewartson II (DS II) equation describes the complex amplitude of surface waves in 



cn ' shallow water 

iut = —'^dlu — 2dzU — dzd~'^\u\'^ — dz^dz\u\'^, 

and was proved in ^23j to be completely integrable. The modified NV equation (mNV) is a 
/\ , member of the DS II hierarchy: 

-■■ Ut = -d'zU-d\-^idzU)(dzd;'\u\')-^(dzU)(dzd:'\u\') 

3 — — 3 — — 

--udzd~^{udzu) - -ud~^{dz{udzu). 

Here dz and d^^ are the solid Cauchy transforms defined by 

d~^u = - -dC, d;^u = - - — =dC 

The integrability of the NV equation has been recently proved in [21] where it is shown 
that a Miura-type map takes solutions of the mNV equation to solutions of the NV equation 
with initial data of conductivity type. This type of initial condition for the NV equation 
was first studied in [15] where it was shown that the inverse scattering method for the NV 
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equation is well-posed for initial conditions of conductivity type. In [T6J it was shown that 
an initially radially-symmetric conductivity-type potential evolved under the ISM does not 
have exceptional points and is itself of conductivity-type. In [T7] evolutions of rotationally 
symmetric, compactly supported initial data of conductivity type computed from a numerical 
implementation of the inverse scattering method for NV are compared to evolutions of the NV 
computed from a semi-implicit finite-difference discretization of NV and are found to agree 
with high precision. This supported the integrability conjecture that was then established 
in [21] where the class of initial data was enlarged by applying Miura-map techniques rather 
than the scattering maps studied in [151 [16], [1^ . In [19j it is shown that the set of conductivity 
type potentials is unstable under C^ perturbations. 

If we consider real solutions u{x,y,t) and let f{x,y,t) = v{x,y,t) + iw{x,y,t), the NV 
equation has an equivalent representation in (x, ?/)-space: 

AUt = -Uxxx + ^Uxyy + 3{uv)x + 3{uw)y, (1.1) 

Ux = Vx- Wy, (1.2) 

Uy = -W^ - Vy. (1.3) 

If the functions u, v and w are not dependent on y, the NV equation reduces to a KdV-type 
equation 

= iUt + Uxxx — GuUx, u = v, w = 

and admits soliton solutions of the form 

u{x,y,t) = — 2c sech.'^ {y/c{x — ct)), (1.4) 

v{x,y,t) = — 2c sech^ {^/c{x — ct))^ (1.5) 

w{x,y,t) = 0. (1.6) 

To facilitate the investigation of the qualitative nature of solutions to the NV equation, we 
present a version of a semi-implicit pseudo-spectral numerical scheme introduced by Feng 
et al [7] that solves the Cauchy problem for the NV equation. This constitutes the first 
numerical implementation of a spectral method for a system of soliton nonlinear PDE's. 
It is shown to preserve the L2 norm for the KdV type equation and is considerably faster 
and requires less computer memory allocation than the finite difference scheme introduced 
in [T7] for the solution of the Cauchy problem. In [1] the method in [7] was applied to the 
KP equation to investigate the stability of soliton solutions. Here, we adapt the method 
for systems of equations to study numerically the nature of the instability of traveling wave 
solutions to the NV equation to transverse perturbations. 

The problem of transverse stability of traveling wave solutions has been studied for many 
of the classic soliton equations including the KP equation [21 El El 113] , the Boussinesq equa- 
tion [5], the ZK equation [21 [HI HDl [H], and most notably, the KdV equation [14j. For the 
NV equation, we carry out a linear stability analysis by considering sinusoidal perturbations 
with wavefront perpendicular to the direction of propagation. In order to draw conclusions 
about the instability of soliton solutions, as well as approximate the growth rate, we ap- 
ply the method developed by Rowlands, Infeld, and Allen [21 112]. Due to the complicated 
boundary conditions, we employ a geometric optics limit based on a scheme that assumes 
the nonlinear wave undergoes a long-wavelength perturbation. Thus, if the wave vector of 
the perturbation is k, we assume it is very small in comparison to the wave vector of the 
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solution. This type of investigation is known as the K-expansion method. To our knowledge, 
the only other use of this method for a soliton system is in the work of Bradley [S] in a model 
of small amplitude long waves traveling over the surface of thin current-carrying metal film. 

In [2], the authors conjecture that to do a linear stability analysis using the K-expansion 
method, a regular perturbation analysis is consistent if and only if the equation is an inte- 
grable system. If a multiscale analysis is needed, the equation is not integrable. The ZK 
equation and the equation in |5J are not integrable, and the ordinary perturbation analysis 
fails. So far this conjecture is supported by the KP and Boussinesq equations [3]. The 
results here support the conjecture by showing that only an ordinary perturbation analysis 
is needed for the NV equation. 

The paper is organized as follows. In Section 2 we show that planar solutions to the 
NV equation must be a solution of a KdV-type equation. In Section 3 the semi-implicit 
pseudo-spectral method is presented, with its linear stability analysis in Section 3.1. The 
K-expansion method is used to establish the instability of traveling wave solutions of NV 
to transverse perturbations in Section 4. Numerical results are found in Section 5 and 
conclusions in Section 6. 



2. From planar solutions of Novikov-Veselov to KdV 

We examine planar solutions to the Novikov-Veselov equations fll.ip - fll.3l) . i.e. solutions 
that depend only on one spatial variable s = (ni,n2) ■ {x,y), moving in a direction given by 
the vector n = (ni, 712) = (cos(q;), sin(a)). We seek solutions of the form 

u(t,s) = u(t,x,y) = u{t,nis,n2s) 

Vi{t,s) = Vi{t,x,y) = Vi{t,nis,n2s) fori = 1,2 

d 

u'{t,s) = --u{t,s) =niUx{t,x,y) + n2Uy{t,x,y) 

OS 

The assumption that u, vi and f 2 are independent on n2X — riiy is equivalent to 

dvi dvi dv2 <9f2 „ 

n2Ux - riiUy = n2 -^ Ui — — = ^2 -^ rii — — = 

ox Oy ox Oy 

As a consequence we obtain u^ = ni u' and q^ = 1^2 u'. The goal is to find a PDE for u{t, s). 
For a given u{t, s) the d equation 
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The nonlinear expression in f ll.ip leads to 

div(-u v) = {nj - nj) (■u^)^. - 2 ni 722 iu^)y + CiU^ + C2Uy 

= (2 {nl — nl) rii — Aril nl) uu' + (ci rii + C2 ^2) u' 
= 2 Kuu + f3u 

where 

/3 = ci rii + C2 n2 and k = K(a) = {nl — nl) ni —2nin\ = cos(3 a) 

The factor k = cos(3 a) is further evidence for the threefold rotational invariance of the 
Novikov-Veselov equations. Now we examine the Novikov-Veselov equation 

13 3 

ut = --u^^^ + -u^yy + -div{uv) 

"^ /// 2 III I n I 

= —-n,u + -ninr,u + - kuu + - pu 
4 i 4 ^ 4 4 

III I n I 

= KU H KUU -\ pu 

AAA' 
Thus, a planar solution of the Novikov-Veselov has to be a solution of the KdV-like equation 

-Ui = -u"' + Quu' +—u' (2.1) 

The only essential modification is the contribution proportional to u' and thus it should not 
come as a surprise that the solutions are related. If q{t, s) is a solution of the standard KdV 
equation 

q(t, x) = —q"'{t, x) + 6 q{t, x) q'{t, x) (2.2) 

then 

u{t,s) =q{-t,s + —t) 

is a solution of the Novikov-Veselov equation (II. ip . which can be verified as follows: 

K 3/3 , 

4 ^ ,. . , 3/3 , ^3/3,3/3, 
—Ut + u — buu u = qt -\ Q + Q — bqq q 

Ho Hh Kj K 

= (lt + q" -Qqq' = 

For a solution g(t, s) of KdV choose constants ki and ^2 with 

3/3 , , 3k 
— = h + k2 — 



and verify that 



u{t, s) = q{- 1 , s + kit) - k2 
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is a solution of the Novikov-Veselov equation: 



K 



Ut 



qt + h q' 



4 , /// « , 3/3 , kiA ,„ , X / 3/3 , 
— Ut + u —buu u = qt -\ q +q —0[q — k2)q q 

AC Av AC Kj 

, fh4: , 3/3\ , 
= qt + q"'-6qq'+{^ + 6k2- — )q' 

= qt + q -6qq + - [h + h— T ) ^ ^ ^ 

Remarks: 

• Consider the choice /3 = ki = k2 = 0. Then q(t,s) = u{jt,s) is the KdV solution 
where the time scale is multiplied with k/4. Since k = K{a) = cos(3 a) the speed of 
the solution changes according to an angle dependent profile, as shown in Figure [TJ 




Figure 1 . Speed profile for planar solutions of the Novikov-Veselov equation 

• The additive constant k2 moves the KdV solution up or down in the graph. 

• Replacing s by s + kit (with ki = ^ — ^2^) corresponds to observing the KdV 
solution in a moving frame, where the frame moves with velocity —ki. 

3. A Pseudo-Spectral method for the solution of (2+1) nonlinear wave 

EQUATIONS 



To numerically solve (ll.ip - (ll.3p we use a semi-implicit leap-frog spectral method based 
on the method in [7]. We restrict ourselves to a finite spatial domain fl = [0, W^] x [0, Wy] 
with periodic boundary conditions. Thus we work on a torus topology and for the numerical 
results we have to observe the traveling waves across the boundary. We must also choose 
the domain large enough for the problem at hand. 

This method uses the Fast Fourier Transform (FFT) to compute the spatial evolution and 
a leapfrog scheme to compute the time-stepping. The numerical scheme can be summarized 
as follows: 

The Fourier transform of (ll.ll) - (ll.3p is 

Aut = i {^^ -3^7]'^) u + 3i^J^[{uv)] + Sir] T[{uw)], (3.1) 

^u = iv-rfw, (3.2) 

(3.3) 



■qu 



-rjv 



uw. 
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As usual, u refers to the Fourier transform F oi u and is given by 

L-l M-1 

1=0 m=0 

L/2-1 Af/2-1 

p=-L/2q=-M/2 

The parameters L and M are the number of grid points in [0, Wx] x [0, Wy\, respectively. They 
need to be powers of two in order to use the standard FFT. The spatial grid is defined by 
(xi, Urn) = (/Ax, mAy) where / = 0, . . . , L — 1, and m = 0, . . . ,M — 1. The spectral variables 
are (^p, V,) = (2vrp/W^„ 2nq/Wy), with p = -L/2, . . , -1, 0, 1, . . . , L/2, and q = -M/2, . . , 
-l,0,l,...,M/2. 

Equations (13.21) and (13. 3p can be solved in terms of u 

Special attention has to be paid to the case E,'^ + 'rf = 0. The corresponding Fourier coefficient 
represent the average values of the functions v and w. 

Let c(t) be the vector with the LM Fourier coefficients of the solution u{x,y,t). Then 
solving (13. 2p . (13. 3p and computing —3i {^J^[{uv)] + riJ^[{uw)]) may be written as one non- 
linear function F{c). With the well constructed diagonal matrix D the NV equation (13. ip 
reads as 

4|c(t)=Dc(t) + F(c(t)) (3.5) 

For the time integration we use a symmetric three-level difference method for the linear 
terms, and a leapfrog method for the nonlinear terms. For the parameter < ^ < 1 we use 
a superscript to refer to the time iteration and the above leads to 

A 

n+l ;?'^— In OT\ /';?'^+l I ;?'^~1 



[C — C 



e D (c"+^ + c"-^) + (1 - 2 0) D c"^ + Fie"") 



2 At 
(2-^AtD)c"+i = (l-20)AtDc'^ + (2 + ^AtD)c"-^ (3.6) 

+AtF(c") 

Thus we have an implicit scheme for the linear contribution and an explicit scheme for the 
nonlinear contribution. Since the matrix D is diagonal we do not have to solve a system of 
linear equations at each time step. For the special case 6' = 1/2 we obtain a Crank Nicolson 
scheme 

(4 - At D) c"+^ = +(4 + At D) c"'^ + 2 At F(c") 

For the three level method we need a separate method for the first time step. For the sake 
of simplicity we may choose c~^ = c°. For computations with known solutions, e.g. (I1.4p . 
we use the known values of the solution at time —At. 



TRANSVERSE INSTABILITY OF SOLUTIONS OF THE NV EQUATION 



3.1. Linear Numerical Stability Analysis. To gain insight into the stability of the spec- 
tral method, we include a linear stability analysis. The parameters determined here were 
used in the numerical experiments that follow. We examine the problem on a domain 
fi = [0, W] X [0, W] C M^ and use periodic boundary conditions. Thus we examine the 
initial boundary value problem 



Aut 

Ur,. = 



Uxxx + 3 Ua:w + 3 a (Vx + w. 



Vr. 



w. 



yi 



"xyy 



yii 



-Wx 



(3.7) 



The case a = is the linearization of f ll.ip - fll.3l) about the zero solution. The utility of a 
linear stability analysis for a nonlinear system was addressed in [7] in the context of the KP 
and ZK equations, where they found that their results for the linear stability analysis were 
validated by numerical results and argued that while the analysis does not prove stability 
and convergence of the nonlinear scheme, the obtained stability conditions often suffice in 
practice. When implemented on a domain with periodic boundary condition (13. 7p leads to 



Uxxx + 3 Uxyy + 3a{vx + Wy), {x,y) xt eQ X 



4ut-- 

Ux = Vx-Wy, Uy= -Wx 

u{x,y,0) = uo{x,y), 

u{x, y, t) = u{x + Wx, y, t), 

u{x,y,t) =u{x,y + Wy,t), 



Using Fourier series and (13. 4p this leads to 

. Sua ^ 
= lUJU + i^^ —u 



e 



v' + C- 



{x,y) e Q, 
{x,y) X te 
{x,y) X te 



u + Sirja 



(3.8) 



-2r7e 



v'+e 



■u 



where u = C,^ — ^r]"^ 



^2 + ^2 

This is an ODE of the form 



4 — u{t) = i A u{t) + i 7 u{t) 

(JjV 



(3.9) 



where A 



OJ 



^^ — S^r]'^ and 7 






Using the numerical scheme (13. 6 p this leads to 



{2-ieXAt) m"+^ = i (1 - 2 0) a At u" + (2 + z ^ a At) n""^ + z 7 At m" 

or with 

61 = 2 + i ^ A At G C, 
to the iteration matrix 

,n+l 



b2=i{(l-29)\ + j)AtEC 



(3.10) 



M" 



U" 



M 



u 



u'- 

n-1 



r b2 


bi -1 


bi 


b^ 


1 


J 



u 



u'- 

n-1 



For the system to be stable we have to verify that the norm of the eigenvalues Zi^2 of M are 
less or equal to 1. The characteristic equation is given by 



^NV{Z) 



— 2- — = 



Thus we conclude |2;j ■ 22! = |t^| = 1 and stability of the solution is equivalent to \zi^2 



R. CROKE, J. L. MUELLER AND A. STAHEL 
61 

i.e. we have |m"'"''^| = lu""^!. Using the explicit solution for quadratic equations and 

b2EiR, |62|<2|6i| (3.11) 

we find 



261^1,2 = 62±V^2 + 4&i&i = &2±-\/&i + 4|6i|2 

l^f'll 1 2^1,2 1 — p2| ^ |02| +4|0i| — \AOi\ 

and consequently 1^1,2 1 = 1 . With f l3Jn|) the stability condition fl3TT|l for the ODE ([33]) is 

((l-2^)A + 7)2(At)2 < 4|2 + z^AAtp 
(((1-2 0) A + 7)2 -(2^ A)') (At)' < 16 
((l-40)A' + 2(l-20)A7 + 72)(At)2 < 16 

For 7 = a = and \ < < 1 the stability condition is satisfied, independent on the step 
size < At and we have a stability result for the initial boundary value problem 



Theorem 1. Apply the numerical scheme (13.61) to the linearization of (ll.ip - (ll.3p about the 
zero solution, i.e. to the equation 

4 xij "Uxxx ' "J "'xyy 

Ux = Vx - Wy , Uy = -Wx - Vy 

This scheme is .stable for \< <1. 

Proving stability for (13. 8 p with a 7^ requires a few more computations. Use the ex- 
pressions for A = a; = ^'^ — 3 ,^ 77' and 7 = 0; -^fn in (13. 9p and the stability condition reads 
as 

((l-40)A2 + 2(l-20)A7 + 72)(At)2 < 16 

3 0- q o-^ 

a;M(l-4 0)+2(l-2 0)-^^^ + ^-^^^)(At)2 < 16 

3a /3a\ /'4^ 



;i-40) +2(1-20)— —+—— < 



^2 + ^2 y^2^^2j - y^/\t 



For = i we find the sufficient conditions. 



3a A / 3a A 

< —. ^ < 



v^+cj ~ w+^y ~ V(^^-3^r72)At 



_3a(^3_3^^2j 

Since |^^ — 3^77^1 < [Sf^ + rf) ■ max{|^|, |?7|} we have the sufficient condition 

4 

At < ^ ,,„ , ,, (3.12) 

3|a| maxjl^l, |?7|| 
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For the DFT with Ax = Ay = ^ we find niax{|^|, |^|} = ;^ and thus we have the stabihty 

condition 

4 
At < ^-^ Ax (3.13) 

As a consequence we have a stabihty result for the initial boundary value problem f l3.8p with 
Theorem 2. Apply the numerical scheme (13. 6p to (13.81) . i.e. 

AUt = -U^xx + 3 Uxyy + 2,a{Vx + Wy) 
UX = Vx-Wy , Uy = -Wx - Vy 

This schem,e is stable for 6 = -^ if the tiwie step At satisfies condition (13.131) . 

4. Instability of traveling-wave solutions of the NV-equation to 

transverse perturbations 

The K-expansion method presented by Allen and Rowlands considers long wavelength 
perturbations in the transversal direction. It was originally used to investigate the stability 
of solutions to the Zakharov-Kuznetsov equation, a non-integrable generalization of the KdV 
equation [TOl [HI |9|. Since then the method has been applied to the KP equation and the 
modified ZK (mZK) equation [18] . This is the second application of the method to a system 
of PDE's, the first being by Bradley [S] . 

4.1. The modified equations. To begin, we transform the system to move along with the 
soliton by using new independent and dependent variables. Using 

{t,x,y) ^ y/c{at,x-ct,y) = {i,x,y) 

and 

{u,Vi,V2) ^ c{u,Vi,V2) 

or 

u{t,x,y) = cu{tc^/^,{x-ct)c^^'^,yc^/^) 

vi{t,x,y) = cvi{t^'\{x-ct)^'\y^'^) (4.1) 

V2{t,x,y) = cv2{tc^/Mx-ct)c'/\yc'/^), 
the NV equations (ll.ll) - (ll.3p are transformed to 

= 4:Ut-4:Ux+Uxxx-3Uxyy-3{uv)x-3{uw)y, (4.2) 

Ux = +Vx - Wy, (4.3) 

Uy = -Wx - Vy (4.4) 

For sake of a more readable notation we dropped the tildes on the new dependent and 
independent variables. The known solution 

fc{x — ct) = —2c sech^(A/c (x — ct)) 

of the original system turns into a stationary solution of the modified system. 

■^0(3;) = t'o(x) = —2 sech^(x), ^0(2^) = 
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Now we examine perturbed functions of the form 

u{x,y,t) = uo{x)+ef{x)e'^y^^\ (4.5) 

v{x,y,t) = VQ{x) + eg{x)e'^y+^\ (4.6) 

w{x,y,t) = wo{x)+eh{x)e'''y+"'\ (4.7) 

Thus we examine periodic, transversal perturbations with exponential growth. For numerical 
purposes we may also work with a purely real formulation. 

u{x,y,t) = uo{x) + ef{x) cos{ky)e'\ (4.8) 

v{x,y,t) = vo{x) + eg{x) cos{ky) e"^^, (4.9) 

w{x,y,t) = wo{x) + e h{x) sm{k y) e^\ (4-10) 

Substituting these expressions into (I4.2p -( B^ and dropping terms propotional to e^ leads 
to (use ' for ^) 

= 47/-4/' + r + 3A;V 

-3{uogy-3{vofy-3uotkh (4.11) 

/' = +g'-ikh (4.12) 

ikf = -h'-ikg (4.13) 

This is a system of linear ordinary differential equations with non-constant coefficients. Use 
h{x) = j^ {g'{x) — f'{x)) and the above can be written in the form 

fix) = (-47 + 3n'o(x))/(x) + (4-3fc2)/'(x) + 

+3uQ{x)g{x) + Quo{x)g'{x) (4.14) 

g"{x) = eg{x) + f"{x) + ef{x) (4.15) 

4.2. Known solutions. Since the fc{x — ct + s) is a solution of the original NV equa- 
tion (ll.ip - 01.3p . we conclude that {uo{x — s), Vo{x — s), Wo{x — s)) solve the system (14. 2p - 
(14. 4p and thus their derivatives must solve the linearized problem (I4.1ip -f l4.13p . Since 
^ sech^(x) = —2 sech^(a;) tanh(a;) we find solutions 

f{x) = g{x) = sech^(a;) tanh(a;) , h{x) = (4.16) 

for /c = 7 = to the system (I4.14p -f l4.15p . For k = 1 and 7 = we found 

f{x) = -|-sech^(a;) 

g{x) = — sech(a;) tanh^(x) (4-17) 

h{x) = —i sech(a;) tanh(a;) 

as a second, explicit solution of (l4.1ip -( H:.13p . 

With the help of these solutions we will construct unstable solution to the Novikov-Veselov 
equations, as function of the parameters k and 7. Thus we need nonzero solutions of (I4.14p - 
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4.3. Locating unstable solutions. Using the function Uq{x) = —2 sech^(a;), a matrix 
notation translates the high order system fl4.14[) - fl4.15[) into a system of first order ordinary 
differential equations. 



d 
dx 



I f \ 




f 




f" 
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9 




\3' ) 








1 

















1 








47 + 3m[) 


(4-3A;2) 





3m'o(x) 


6 Mo 
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fc2 





1 


fc2 






/ / \ 

/' 
/" 

9 
\9' J 



4.3.1. Behavior of solution for large \x\. For large values of |a;| we use uo{x) ~ and arrive 
at a decoupled system of linear differential equations with constant coefficients. 



(4.18) 
(4.19) 




/' = 

f" J 





-47 


1 
1 

(4-3A;2) 


f 
\f" 


3\ 
9' J ~ 


" 1 " 
k^ 


f 9\ ( 

, + 
\9 V 




ef' + 



The characteristic equation of the matrix in (I4.18P is given by 



(4.20) 



and we will denote the eigenvalues Ai, A2, and A3 of this system by pi,p2, andps, respectively. 
The inhomogeneous system in 04.19P leads to eigenvalues A4^5 = ±fc. Using Cardano's 
formulas for zeros of polynomials of degree 3, we verify that 04.2UJ) has three distinct real 
values if 

, 2 (3A;2-4)3 

For our domain of (/c, 7) to be examined we may use 

P3 < < pi, P2 
and for small values of I7I we find 



A(7) 



+ i4li7 

+ v /4^3F -^4^7 

-v/4^3F-^4^7 



(4.21) 



with the solution of the form f{x) ~ cie^^^ + 026^^^ + 036^^^. We examine 7 > 0, and for 
X — )■ +00 we require the solution /(x) of (I4.18P to be bounded, and thus 



f{x) K. /+(x) = cse 



P^x 



Next, examine (I4.19p in the form 



gl{x) - k' g+ix) = k' f+{x) + f'lix) = C3 {k' + pj) e^^' 



with the solution 



9+{x) 



kx 



/3ie+'=" + /32e-'^" + C3 



k"^ +pI 
pI — k"^ 



^P3X 
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Since this solution has to remain bounded as x — )■ +00, we find 

kx , „ ^ +P3 ns X 



g+{x) =(326 '=^• + 03 



pI — k? 



The case 7 > and x — > — 00, leads in an analogous way to 
and 



3l 6+^^^ 


Pi 




ePi" + C2 


fc2 
Pi 


+ pI 


P?-7 


42 

A J 


T^7^ 


^7^ 







pP2a; 



For < A;^ < 1 we find 

/l2 

(4.22) 

and thus the above formula for g(^x) may not be valid for 7^ ^ A;^ <^ 1 . For /c = the equation 
for (yf simplifies to g"{x) = f"{x), and the conditions at a; = ±cx3 imply g{x) = f{x). 



4.3.2. Solutions for intermediate \x\. Choosing a large value for M > 0, we construct non- 
trivial solutions to f l4.14p -f l4.15p for —00 < x < — M, then for —M < x < +M, and then for 
+M < X < -l-cxD. We seek nonzero values of the parameters (ci, C2, C3, /3i, /32), such that we 
find a nonzero solution. 

Use the solutions /_ and g- to define a matrix 



g-pi Af 

p2g-piM 
Pi+fc^ „-pi M 



g-P2 A/ 

p2g-p2Af 

PJ+fe^ P-P2 M 

(p|+fc^)p2 -p2 M 







^-kM 

A;e-"^ 



with 



C2 
C3 

Use the solutions /+ and g^ to define a matrix 

pPsM 



( fi-M) \ 

n-M) 

g{-M) 

\ g'i-M) J 











Pa 



p2gP3M 



^psM 


















-kM 



(P3+L)P3gP3A^ 



-ke 



-kM 
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with 



Pi 



( fi+M) \ 

n+M) 
n+M) 

g{+M) 
V 9'{+M) I 



Then use an ODE solver to examine the system fl4.14p - fl4.15p on the interval [— M, M] as an 
initial value problem. This leads to a matrix T such that 



n-M) 

g{-M) 

\ g'i-M) I 



( fi+M) \ 

n+M) 
n+M) 

g{+M) 
V 9'i+M) I 



Now we have two methods to compute the values of the solution at x 
system fl4.14p - fl4.15p has a nonzero solution if and only if 



+M, and the 



D{k, 7) = det(M(fc, 7)) = det(T ■ T_ - T+) = (4.23) 

Thus we examine solutions of this equation as functions of the parameters k and 7. 

4.4. The special case fc = 0. For fc = 0, equation fl4.15p reads as g"{x) = f"{x) and the 
only solution satisfying g{±oo) = /(±oo) is g{x) = f{x). Then (14.141) leads to 




1 

1 

-47 + 6m'o 4 + 6mo 




(4.24) 



Using the same approach as in the previous section we define 




fi-M) 

n-M) 
n-M) 

fi+M) 

n+M) 
n+M) 

fi-M) 

n-M) 
n-M) 



-piM 



pi e Pi *^ p2 e 



-P2 M 
-P2 M 



pj e-Pi *^ 



p2 g-p2 M 






















gPsM 
p2gP3M 




fi+M) 

n+M) 
n+M) 



where Tq is constructed using an ODE solver for fl4.24p . Then the system (14. 140 - 04. 15p has 
a nonzero solution for k = ii and only if 



Do (7) = det(Mo(7)) = det(To ■ T_ 







(4.25) 
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5. Numerical Results on the Instabilities of Plane- Wave Soliton Solutions 



5.1. Locating unstable solutions. Based on expression fl4.23p we generate plots of the 
function D{k,'~f) = det(M(A;,7)) on a domain < fc < 1 and < 7 < 0.5, leading to 
Figure |2j In the corner k ^ 1 and 7 ~ 0.5 the real part of the function vanishes, but a 
second plot verifies that the imaginary part is different from zero. Thus Figure [2] indicates 
that we have a clearly defined solution curve of det(Wl{k,'y)) = 0, away from the origin. 




0^0 



Figure 2. The real part of det(M) as a function of k and 7 



Since the behavior close to (^,7) ~ (0, 0) is critical, we examine this section with a 

The obvious spikes are caused by the zeros in the 



finer resolution, leading to Figure 3 (a 



denominator in condition fl4.22p . Figure 3(a) suggests the existence of a solution along the 
axis k = 0. Using (14.251) we generate Figure 3(b), As a consequence the only solution along 
the axis fc = is at 7 = . Thus the known solution (14.161) is an isolated solution in the 
parameter space (^,7) at (0, 0). 

With the above preparation we can now construct values of [k, 7) leading to nonzero 
solutions of (I4.14p - (l4.15p and thus for 7 > to unstable soliton solutions of the NV equa- 
tions flTTTIl -f lT:^ . 

We trace a solution curve of (I4.23P by generating an arc length parametrization of the 
curve. Using an arbitrary initial point on the curve (use a contour plot of Figure [2]) and a 
starting direction, we minimize | det(M(/c,7))| along a straight line segment orthogonal to 
the stepping direction. Using this local minimum, we adjust the stepping direction and then 
make a small step to follow the solution curve. While stepping along the curve we verify 
that we actually have a solution of (14.230 . and not only a minimum. This proved to be a 
stable algorithm and generated Figure HI 
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(a) det M(A;,7) close to origin 
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det Mo (7) 





Figure 3. Behavior at the origin and along the axis k = 



Observe that Figure 4(a) also displays negative values of 7, which do not lead to unstable 



solutions of f ll.ip -f lO]) . We display these values to confirm that we have a closed curve 
without branching points. 



O.S r- 
0.4 - 
0.3 - 
0.2 - 
f- 0.1 

-0.1 
-0.2 
-0.3 — 




/ 



0.2 0.4 O.G O.S 



(a) solution curve of | det(M)| — 



(b) parameters leading to unstable solu- 
tions 



Figure 4. Solution curve of |det(M)| 
to unstable solutions 



and parameters k and 7 leading 



For the parameter values of (fc, 7) in Figure 4(b) there are nonzero functions /(x), g{x) and 
h{x), such that for small e we have solutions of the equations (14. 2p - (14. 41) of the form (14. 5p - 
(14 .yp . Since the y dependence of these functions is of the form e^^^, the functions are periodic 
in y with a period oi Ly = ^. The corresponding exponent is shown in Figure (4(b) ). Thus 
the Novikov-Veselov equations (ll.ll) - (ll.3p with initial condition u{x, y, 0) = — 2 c sech.'^{y/cx) 
will lead to an unstable soliton solution. 

Theorem 3. The soliton solutions u{x,y,t) = v{x,y,t) = — 2csech (a/c(x — ct)) and 
w{x,y,t) = of the NV equations (ll.ip - (ll.3p are not stable. For values 0/ 0.363 < A; < 1 
there are y -periodic, unstable contributions with a period of Ly = y%- 
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5.2. Constructing unstable solutions numerically. As an example, in this section we 
construct one of the above unstable solutions numerically, using the algorithm from Section [3] 
for periodic solutions in x and y. Since the soliton solution sech^(x) decays rapidly, the above 
results still apply when working on a sufficiently large domain. The steps of the algorithm 
are as follows: 



(1) Choose values of {k, 7) in Figure 4(b) and determine the eigenvector (ci, C2, C3, /3i, /32) 
for the zero eigenvalue. 

(2) Use the algorithm leading to the matrices T, T_ and T+ to construct the nonzero 
functions f (x) , g (x) , h{x) . 

(3) Pick a size domain such that ^ periodic functions in y are admissable. 



(4) Construct initial values, using f l4.8p - fl4.10p . use a small value of e. Without taking 
the transformations (14. ip into account we obtain a speed of c = 1 of the unperturbed 
soliton. 

(5) Solve the NV equations fll.ip -f lL3|) . using the algorithm presented in Section |3l 

(6) The deviation from the single soliton solution should not change its shape, but the 
size is expected to be proportional to e'*'*. 

(7) Sohtions for speeds 07^! can be constructed similarly, using the transformations f l4.ip . 

The evolution of a perturbed soliton from Theorem [3] was computed using the semi-implicit 
pseudo-spectal method. Here we chose k = 0.504 and 7 = 0.296, and the graph of the 
perturbation f{x) is found in Figure 6(a) As initial value we chose a perturbed KdV 
solition with speed c = 1, starting at x = 10. Find the solution and the difference to 
the unperturbed KdV soliton at time 5 in Figure |5l The corresponding animations are 
available on the web site [22] . The exponential growth of the perturbation with exponent 7 
is numerically confirmed. 



Solution 



Difference 













(a) shape of the solution 



(b) shape of the perturbation 



Figure 5. A perturbed solution at time 5 



One can construct the shapes of the functions f{x) for all positive values of 7 along the 
arc in Figure H] to obtain Figure 6(b)[ The solutions constructed for (A;,7) = (1, 0) have to 
match the known exact solutions (14.170 . which is confirmed. 
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0.4 

0.2 



-0.2 

-0.4 




(a) at fc ~ 0.5 and 7 w 0.3 



(b) as function of arc length s 



Figure 6. Shape of the perturbations f{x) 



6. Conclusions 

In this work a semi-imphcit pseudo-spectral method was introduced for the numerical 
computation of evolutions of solutions to the NV equation, constituting the first numerical 
implementation of a spectral method for a system of soliton nonlinear PDE's. A linear 
stability analysis yields a stability condition for the Crank-Nicolson scheme on the linearized 
IBVP. The instability of traveling wave solutions to transverse perturbations was established 
by the K-expansion method, and unstable soliton solutions were constructed. The evolution 
of an example was computed numerically by the semi-implicit pseudo-spectral method. 
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